function[y] = F(p_star)
global gamma p1 rho1 c1 p2 rho2 c2 u1 u2;

if(p_star >= p1)
    f1 = (p_star - p1)/rho1/c1/( (gamma+1)/2/gamma *(p_star/p1) + (gamma-1)/2/gamma)^0.5;
end
if(p_star <p1)
    f1 = 2*c1/(gamma-1)*((p_star/p1)^((gamma-1)/2/gamma) - 1);
end
if(p_star>=p2)
    f2 = (p_star - p2)/rho2/c2/( (gamma+1)/2/gamma *(p_star/p2) + (gamma-1)/2/gamma)^0.5;
end
if(p_star <p2)
    f2 = 2*c2/(gamma-1)*((p_star/p2)^((gamma-1)/2/gamma) - 1);
end
y = f1+f2 + u2 -u1;
end